13.3 s
2.5 ms
2.7 μs
2.4 ms

Modeling and Optimization

An epidemic example

6.5 μs

Want all the details?

This talk is based on the free, online course "Computational Thinking with Julia" available from juliaacademy.com


Lecture notes and problem sets (.ipynb and .pdf) are available at https://github.com/mitmath/6S083/


An updated and expanded version of the course called "Introduction to Computational Thinking" is available at https://computationalthinking.mit.edu/Spring21/

15.2 μs

83.0 ns

Further reading:

2.7 μs

67.0 ns

Algorithms for Optimization by Kochenderfer & Wheeler

5.4 μs

73.0 ns

Algorithms for Decision Making by Kochenderfer, Wheeler, & Wray

4.9 μs

The Procedure

1. Develop a computer model to simulate reality

2. Establish a theoretical foundation to explain results

3. Compare the computer model to the theory

4.9 μs

The Random Walk Model

Let N "agents" randomly move within a square of side L.


For true "agent-based modeling" (ABM), use Agents.jl.

11.6 μs

First define the walkers and their methods.

2.6 μs

Use an "enumerated type" to keep track of infection status

2.3 μs

Susceptible -> Infectious -> Recovered

8.1 μs
7.1 ms
R::InfectionStatus = 2
76.0 ns
173 μs
52.5 μs
40.6 μs
65.6 μs
walkers
3.8 ms





32.0 ns





79.0 ns
28.7 ms
20000
23.2 ms
43.2 μs

The Simulation

2.3 μs

Define close contacts using NearestNeighbors.jl

2.0 μs
63.9 μs





52.0 ns

step! sweep! sim!

4.8 μs

One step! moves one agent in one time step.

One sweep! moves all agents on average in one time step

One simulation is a series of time sweeps

2.6 ms
54.8 μs
26.5 μs
64.5 μs

Visualize it

6.2 μs
48.2 μs
4.0 s
75
30.6 ms
6.3 ms





83.0 ns
117 ns
sim
24.1 s
128 ms

The SIR Model

2.9 μs

Has three states:

  • Susceptible

  • Infectious

  • Recovered

The time evolution of these states can be described by a set of ordinary differential equations.

3.8 ms

dSdt=βS(t)I(t)

48.9 ms

dIdt=βS(t)I(t)γI(t)

32.5 μs

dRdt=γI(t)

8.2 μs

48.0 ns

The approximate solution can be taken in small time steps.

6.0 μs

S(t+h)S(t)hβS(t)I(t)

5.7 μs

I(t+h)I(t)+h[βS(t)I(t)γI(t)]

2.7 μs

R(t+h)(t)+hγI(t)

6.3 μs
279 μs
test
43.6 μs
18.7 ms

Optimizing with Gradient Descent

6.3 μs
1.4 s

Use derivatives of a function to find its mininmum.

5.8 μs





82.0 ns

The derivative of a function can be approximated as:

6.0 μs

dfdx(a)f(a+h)f(a)h

6.2 μs
f (generic function with 1 method)
23.5 μs
deriv (generic function with 2 methods)
29.5 μs
tangent_line (generic function with 1 method)
165 μs
2.1 ms

a = 3

79.2 μs





31.0 ns

Extending to 2-D with the gradient, ∇

3.8 μs

For a function g(x, y), the gradient is a vector of partial derivatives with respect to x and y.

4.7 μs

g(x,y)=x2+y3

6.0 μs

gx=2x

2.8 μs

gy=3y

3.4 μs

g(x,y)=(g/xg/y)=(2x3y)

4.4 μs

g(2,5)=(2(2)3(5))=(475)

4.6 μs

34.0 ns
20.7 μs
61.4 μs
59.2 μs
48.6 μs
9.9 ms





73.0 ns

Gradient descent optimization in 1-D

3.4 μs
135 μs





75.0 ns

Let's try it out!

3.0 μs
22.7 μs
x_descent
383 ms
y_descent
24.3 μs
2.8 ms

x₀ =

89.7 ms

index = 1 326

7.8 ms





87.0 ns

Carl Friedrich Gauss walks into a bar...

3.4 μs

and tries to fit his data.

3.2 μs

33.0 ns

G=Aσ2πe(xμ)22σ2

6.1 μs
33.8 μs
27.7 μs
992 ms





69.0 ns

Gradient descent optimization in 2-D

2.5 μs
78.2 μs

Use gradient descent to "learn" μ and σ

2.4 μs
82.0 ns
xg
-10.0:0.005:10.0
2.5 μs
yg
71.7 ms
264 ms
125 ms





82.0 ns

The loss function measures how much the model deviates from the data.

3.7 μs
37.8 μs
2.5 μs
2.1 μs
z (generic function with 1 method)
19.4 μs

34.0 ns
902 ms





124 ns

Passing loss into gradient_descent selects parameters that drive it towards zero.

3.8 μs
959 ms
32.1 ms

step = 90000

14.6 ms

How does our SIR model measure up?

2.6 μs

Need to redo the loss function

2.6 μs
46.7 μs
30.4 μs
ode_params
3186×2 Matrix{Float64}:
 0.1     0.0
 0.1001  0.0001
 0.1002  0.0002
 0.1003  0.0003
 0.1004  0.0004
 0.1005  0.0005
 0.1006  0.0006
 ⋮       
 0.418   0.0164
 0.4181  0.0163
 0.4182  0.0164
 0.4183  0.0163
 0.4184  0.0164
 0.4185  0.0163
464 ms
86.2 μs





86.0 ns
4.0 ms

Fit Parameters

β = 0.4185

γ = 0.0163

49.3 μs

step = 1 3186

54.5 μs

β = γ =

17.6 ms